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Abstract 

First-principles atomistic molecular dynamics simulation in the micro-canonical and canonical 
ensembles has been used to study the diffusion of interstitial hydrogen in a-iron. Hydrogen to 
Iron ratios between = and \ have been considered by locating interstitial hydrogen atoms 
at random positions in a 2 x 2 x 2 supercell. We find that the average optimum absorption site 
and the barrier for diffusion depend on the concentration of interestitials. Iron Debye temperature 
decreases monotonically for increasing concentration of interstitial hydrogen, proving that iron- 
iron interatomic potential is significantly weakened in the presence of a large number of diffusing 
hydrogen atoms. 
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Introduction. Diffusion of hydrogen atoms in metals make an impact on a number of 
technologically important properties of the host material and has attracted interest since 
long ago.- In particular, it is believed to play a crucial role in embrittlement of high strength 
steels, of common use in buildings, bridges, etc. 2-5 However, a fundamental understanding of 
different mechanisms linking embrittlement and diffusion of hydrogen is still lacking owing 
to the complexity of the relevant processes. Different authors have used Density Functional 
Theory (DFT) to study preferred adsorption sites and diffusion barriers for hydrogen in 
body-center cubic iron (a-Fe), which constitutes the core of high-strength steels and con- 
sequently is the simplest model for fundamental studies trying to connect macroscopic and 
microscopic properties.- - - The main stream in the literature favors diffusion barriers of ~ 0.1 
eV, and adsorption of hydrogen in the tetrahedral site which is the mot favorable one simply 
in terms of available space. Comprehensive work along these lines, including the effect of 
surfaces, has been done by Jiang and Carter, and and more recently Ramasubramaniam 
et al.— ^ For large concentrations, however, we have found that occupancy of the compet- 
ing octahedral site, and reduced diffusion barriers, can happen.— This is accompanied by 
a tetragonal distortion in the bcc unit cell and it is related to internal stresses appearing 
upon an increasing hydrogen load inside the system. From an experimental point of view, 
it has been observed that for increasing external pressures of hydrogen the iron melting 
point is lowered appreciably and the transition from a (bcc) to 7 (fee) happens at lower 
temperature.— Furthermore, neutron diffraction analysis on ferrite samples of high-strenght 
steels have revealed a significant increase of Debye- Waller factors upon increasing loads of 
hydrogen. 14 In this paper we investigate the effect of accumulating interstitial atoms in a 
region, which is known to happen at least locally around defects prior to the breaking of 
samples. We use ab-initio molecular dynamics because it presents several advantages to 
study this problem. First, it allows us to study the many-body interactions between simul- 
taneously diffusing interstitial impurities. Second, temperature dependent simulations allow 
us to extract useful additional dynamical information. Third, for such a complex situation 
where the number of possibilities to distribute a number of interstitial impurities among sev- 
eral candidate sites becomes combinatorially large, a statistical approach where the system 
follows its own internal dynamics to sample relevant configurations is the more effective ap- 
proach. Kinetic Monte-Carlo has also been recently applied to study the variation of barriers 
with different configurations and their associated stresses:— this is an alternative approach 
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where good agreement with prior ab-initio DFT calculations has been found, proving the 
viability of statistical methods to study this kind of problem. 

Methods. First-principles MD calculations have been performed with the CASTEP 
code.— ^ A 2 x 2 x 2 periodic supercell is set up including sixteen Fe atoms and several H 
atoms depending on the different concentrations being considered. The Born-Oppenheimer 
approximation is used; ions are considered classical objects moving on the potential created 
by the electrons obeying the Schroedinger equation. Electronic wave-functions are expanded 
using a plane-waves basis set up to an energy cutoff of 375 eV and are sampled inside the 
Brillouin zone in a 4 x 4 x 4 Monkhorst-Pack mesh.— Electronic energies are converged up 
to 10~ 5 eV. Ultra-soft pseudopotentials are used to describe Fe and H— and the generalized 
gradients approximation for exchange and correlation due to Perdew, Burke and Ernzerhof 
is chosen^. These approximation have been thoroughly checked before and it has been 
found that they reproduce correctly the main physical properties of a-iron, including lattice 
constant, magnetization and bulk modulus.— A word of caution is in order here since we 
have used a k-point mesh less dense than the ones we have previously shown to be adequate 
to accurately reproduce different equilibrium properties of a-iron at T = K. Our choice 
is based in two different reasons. First of all, it is a practical one since molecular dynamics 
is a computer intensive task and a compromise between accuracy and time must be made. 
Secondly, we notice that computing different physical magnitudes require different accura- 
cies, as can be deduced from Fig. 2 in ref.—; computing the equilibrium lattice parameter 
with a k-point mesh similar to the one used here incurs in an error of 0.4%, but to ~ 300% 
for the bulk modulus. Our primary goal in this work is to compute total energies that can 
fluctuate in the current range of temperatures by ~ 0.05 eV. We have checked that for the 
maximum density of interstitial impurities considered here the total energy changes at T = 
K by 0.047 eV if the 4x4x4 mesh is replaced by a 8 x 8 x 8 one. Therefore, the error is 
close to the random fluctuations intrinsic to the system and can be accepted. Furthermore, 
we notice that the error in the equilibrium lattice parameter by such an approximation is 
rj 0.01 A (e.g. Fig. 2 in ref.—), which is of the same order of magnitude or lower than the 
root mean squared amplitude displacements of vibrating atoms at the temperatures used 
here. Therefore, we assume that the 4 x 4 x 4 is both practical and accurate enough for 
the purpose of these simulations. Finally, our current MD simulations predict a value for 
the Debye temperature of bcc iron (O = 505 K) in reasonable agreement with the reported 
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experimental value (420 K). 

To study the diffusion of several H atoms inside the unit cell we test in the micro-canonical 
ensemble the quality of the total energy conservation during a typical MD run. Simulations 
for 1 or 2 ps are performed with time steps of 0.5 or 1.0 fs showing that the total energy 
is conserved within a 0.01% error (equilibration taking place during the first 100 to 200 fs 
have been taken out from averages). Keeping fixed the parameters defining the model, we 
switch to the canonical ensemble to reproduce conditions where the Fe and H atoms are in 
equilibrium with a thermal bath kept at a fixed temperature (Nose-Hover prescription has 
been used). All these simulations are performed keeping constant the volume of the unit 
cell and the number of particles. These boundary conditions are important to understand 
the physical model and the solutions obtained. In particular, it is relevant to discuss the 
meaning of keeping the volume of simulation fixed. In a previous work we have investigated 
the volume deformation and atomic displacements necessary to find an equilibrium solution 
with zero forces and stresses in the presence of interstitial hydrogen binding to either the 
octahedral or tetrahedral high symmetry sites.— Here we are interested, for the case of 
an overall low dilution concentration of impurities, in the effect of a high concentration 
of interstitial atoms inside a small region embedded in a matrix of iron that except for 
the large concentration of interstitial hydrogen in a small number of regions mostly keeps 
its original properties. This is consistent with the experimental observation of overall low 
dilution of interstitial H in a-iron, but possibly large concentration in particular regions, 
maybe as a consequence of the existence of defects.- 1 ^ Therefore, it is assumed here that 
the modification of the volume in the region of interest due to the internal pressure created 
by the impurities is effectively controlled by the large amount of unperturbed bulk material 
surrounding the region where the interstitial hydrogen diffuses in a scale of time compatible 
with our simulations (ps). Consequently, we fix the 2x2x2 supercell volume used for 
the simulation to the one corresponding to pure a-iron independently of the interstitial 
concentration inside. In the other limit a completely different model could be considered: 
an overall concentration distributed uniformly and large enough large enough so the volume 
should be adjusted for each density. For this case, the fractional variation of the supercell 
volume would be 4%, 7% and 15$ for n H = 2, 4, 8, respectively (the restriction a = (3 = 7 = 
90° has been used, but a, b, and c have been allowed to change freely to minimize stresses). 
While the (NVE) and (NVT) collectives are appropriate for the first scenario above, for 
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the second one the (NPH) and (NPT) collectives should be used. The (NPH) and (NPT) 
collectives would also be probably more adequate to describe the region near the surface, 
where the lattice of iron breaks its periodicity and might not be strong enough to counteract 
the internal pressure due to impurities. This work has been restricted to the (NVE) and 
(NVT) collectives, leaving the used of other collectives for future research. 

Diffusion coefficients are computed by assuming a random-walk for interstitial 
impurities.— Fig. [TJ shows the different paths for the high density case with eight hydro- 
gen atoms diffusing simultaneously inside the 2x2x2 supercell. As already predicted 
by our static geometry optimization interstitial hydrogen atoms avoid each other and no 
tendency to formation of molecular hydrogen has been observed. We follow the different 
trajectories during the simulation time and compute averages 

< \r(t) -r(0)| 2 >= 6Dt (1) 

to extract the three-dimensional diffusion coefficient D. Barriers for diffusion are estimated 
from an Arrhenius plot D = D e~ B ^ kBT by a least-squares fit. Here, the prefactor D is 
related to a typical vibrational frequency for H in the Fe bcc lattice, and it is assumed to 
be a constant value independent of the number of interstitial atoms in the supercell. This 
assumption is corroborated by our fits within an uncertainty of ±12% (Fig. [2]). 

Results. We simulate MD trajectories for a single H atom diffusing inside the supercell 
to compare with previous results derived from ab-initio DFT geometrical optimization and 
transition state theory. From a least-square fit to the data in the Arrhenius plot we deter- 
mine a barrier for diffusion B\ = 0.145 eV (Fig. [2]) in good agreement with previous values 
obtained under similar conditions (cubic symmetry, and the same interstitial density^). 
This agreement shows that our MD simulations adequately sample the relevant phase space. 
The time evolution of the interstitial atom can be further analyzed to show that under this 
conditions trajectories near tetragonal sites (T) are preferred over octahedral sites (O). This 
conclusion can be made quantitative by computing a characteristic residence time for both 
sites. We assign parts of the diffusing path either to T or O according to proximity to esti- 
mate the likelihood to find the particle, say, near O. Fig. |3] displays the fractional occupation 
of octahedral sites for 1, 2, 4 and 8 hydrogens in the 2x2x2 supercell. These values can 
be understood by comparing with a simple two-level model where the only parameter is the 
energy difference between T and O sites, k B A = E — E T : 21 
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+ T 1 + 2e A / T 

This model reduces all the accesible volume for the interstitial hydrogen diffusing in the 
unit cell to only a set of discrete lattice points (either T or O), but in spite of its crudeness 
it already captures the essentials of the problem. Eq. [2] has been plotted in Fig. |3] for 
&bA = 0.06, 0.07, 0.04 and —0.035 eV, yielding least-square fits to the points extracted from 
the MD simulations for uh = 1,2,4 and 8 interstitial hydrogens, respectively. The dashed 
line represents the asymptotic equilibrium distribution (Ot->oo = §) to be approached from 
above or below depending on the sign of A. These results show that for n# = 8 (9 = ~) 
A has moved from positive to negative and the equilibrium site has been exchanged from 
T to O. The dependence of the parameter A with the interstitial density proves how site 
adsorption energies are affected by the presence of an increasing number of hydrogen atoms 
concentrated in a particular region. This behavior follows the pattern previously predicted 
by ab-initio DFT geometry optimization where the T site is the lowest energy configuration 
for low densities, while for large densities occupancy of the O site favors a body centered 
tetragonal distortion of the lattice and becomes the global minimum. Although for the 7 
phase, experiments reporting a qualitative modification of the system around a hydrogen 
concentration of ~ 0.4 show how the increasing density of interstitial impurities might 
significantly modify the dynamics of these systems (e.g., Fig. 13 in^). This is an observation 
that might help to explain the spreading of values extracted from different experimental 
techniques for diffusion barriers that has previously been linked to partial occupation of 
both sites.— Partial occupation of both sites at a given temperature happens most naturally 
in molecular dynamics simulations, but it is not easy to describe in a standard geometry 
optimization. We remark that the present approach represents a feasible route to investigate 
a regime that otherwise is too difficult: in the 2x2x2 there are 48 different O-sites and 
96 different T-sites, being the number of ways to distribute several interstitial among these 
combinatorially large (~ 10 11 ). Such a huge configurational space can only be addressed 
from a statistical point of view and by letting the system to explore as many relevant cases as 
possible by following its own dynamics. To understand to which extent barriers for diffusion 
are affected by the presence of extra interstitial atoms we analyze the high-density case 9 = \ 
in more detail. From an Arrehnius plot (Fig. [2]) we extract by a least-square fit an effective 
barrier of B 8 = 0.047 eV; significantly lower than the one found for a single interstitial, B\. 



Based in ab-initio DFT geometrical optimization we have previously suggested that an 
important consequence of interstitial hydrogen diffusing in a-Fe is to weaken the Fe-Fe 
interaction.— From a physical point of view this effect is related to several factors: first, 
the presence of interstitial impurities increases the effective Fe-Fe distance and screens their 
interaction; second, the symmetry is distorted, an effect that we have found is important 
to explain the stability of octahedral sites in the high density regime; third, spin-polarized 
DFT calculations reveal that hydrogen contributes an extra spin to the system, but the 
total ferromagnetic moment does not grow accordingly. Current MD simulations should be 
taken as a highly controlled and clean theoretical experiment which results can extend these 
interpretations from T = K to a finite temperature, and from a limited amount of atoms 
sitting in the same supercell, to a larger number meeting together even if for a limited amount 
of time. It is interesting to observe that our MD simulations support a similar interpretation 
to the one derived from ab-initio DFT at T = OK. We compute the mean square amplitudes 
of iron atoms vibrating around their equilibrium positions, < u 2 >, which is related to the 
strength of the potential confining these atoms. For each temperature we obtain the mean 
squared amplitude of vibration by fitting their time-averaged probabilities to a Gaussian 
distribution centered around its equilibrium position. These values are compared with an 
isotropic Debye- Waller model for the mean squared displacement of atoms vibrating at 
temperature T:— 



where G is the Debye temperature and M the mass of the atom. Fig. H] shows how the root 
mean squared amplitude of vibration (u) increases steadily with the number of H atoms 
inside the 2x2x2 supercell, being nearly doubled from tie = 1 to 8. Adopting a Lindemann- 
like criterion we can conclude that increasing the number of interstitial hydrogen drives the 
material closer to a thermodynamic instability that eventually should lead either to a phase 
transition or to the material failure. This idea is more clearly illustrated by using Eq. [3] to 
fit these < "u 2 > to an effective Debye temperature for each density (Fig. [5]). G decreases 
monotonously when the number of interstitial hydrogens increases marking the softening of 
en-iron, and proving that the material is, at a given fixed T, getting closer to its own melting 
point under the internal pressure of dissolved H. This is closely related to the growing number 
of interstitial atoms sitting together in the same unit cell and the inverse situation (growing 
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temperature getting closer to the melting point at a fixed number of interstitial impurities) 
cannot be inferred from our simulations because the Debye Temperature is constant under 
these conditions. 

Conclusions. We have found by direct analysis of our MD trajectories that the diffusing 
barrier for interstitial hydrogen inside a-iron depends on the density of diffusing atoms in 
the near region. By using a simple statistical model we have also analyzed how the energy 
difference between the T and O sites, A, is modified by the presence of other interstitials. 
Finally, by looking at the amplitude of vibration of iron atoms around their equilibrium 
position, and comparing with a simple Debye- Waller model, we conclude that the Fe-Fe 
interaction weakens as the concentration of interstitial hydrogen increases, finding that for 
the largest considered density the effective Debye temperature for iron is already below room 
temperature. 

This work has been financed by the Spanish CICYT (MAT2008-1497), and MEC (CON- 
SOLIDERS CSD2007-41 "NANOSELECT" and "SEDUREC"). 
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FIG. 2: (a) (short dashed line) Diffusion coefficient, D (^-), obtained for a set of averaged diffusion 
MD trajectories for a single H atom in a 2 x 2 x 2 bcc iron supercell. A barrier for diffusion of 
Bi = 0.145 eV is obtained from a least-squares fit. (b) (long dashed line) Same for eight H 
interstitial atoms diffusing in the supercell. A barrier of B$ = 0.047 eV is deduced from the fit. 
Error bars have been estimated from standard errors ). 



11 



0(%) 




FIG. 3: As a percentage over the total simulation time, residence time around octahedral sites, 
Po for: (i) a single interstitial hydrogen (small gray dots; thin gray line corresponds to fcgA = 0.06 
eV in Eq. [2]); (ii) two H (small black dots; thin black line is for fcgA = 0.07 eV); (iii) four H (large 
gray dots; thick gray line is for /cgA = 0.04 eV); and (iv) eight H (large black dots; thick black 
line is for fee A = —0.035 eV). The dashed line separates the two asymptotic regions (A > and 
A < 0). 
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FIG. 4: Iron root mean squared displacement from equilibrium positions u(A) vs. temperature 
(K). The density of interstitial hydrogen for each case is: (i) nn = 1 (small gray dots), (ii) tih = 2 
(small black dots), (iii) nn = 4 (large gray dots), and (iv) n# = 8 (large black dots). Corresponding 
continuos lines coded similarly in size and gray/black are least-square fits to the data using Eq. [3] 
with Debye temperatures = 381,399,360 and 259 K, respectively. 
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FIG. 5: Debye- Waller temperature (K) vs. number of interstitial H diffusing in the unit cell. 
Dashed line is a linear least-squares fit to guide the eye. 
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